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ABSTRACT 


This  work  studies  the  use  of  a  wave  absorbing  control  law  for  vibration 
suppression  of  flexible  spacecraft  structures.  A  major  advantage  of  this  method  is  that  it 
does  not  involve  truncation  into  a  finite  dimensional  mathematical  model  A  closed  loop 
scattering  matrix  was  derived  which  gives  the  relationship  between  incoming  waves, 
outgoing  waves,  sensor,  and  actuator.  The  control  law  was  determined  by  minimizing  the 
H-infinity  norm  of  this  matrix.  The  control  law  was  applied  to  the  Naval  Postgraduate 
School's  Flexible  Spacecraft  Simulator  (FSS)  for  vibration  suppression.  The  simulator's 
flexible  beam  was  controlled  using  piezoelectric  ceramic  wafers  as  sensors  and  actuators 
The  H-infinity  wave  absorbing  controller  contributed  significant  damping  to  the  structure, 
especially  at  the  first  mode  of  1  Hz  Therefore,  wave  absorbing  control  and  piezoceramic 
sensors  and  actuators  offer  a  viable  approach  for  vibration  suppression  of  space  structures. 
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1.  INTRODUCTION 


A.  BACKGROUND 

Recent  trends  in  spacecraft  design  have  led  to  larger  and  more  flexible  structures. 
Large  structures  such  as  Space  Station  Freedom  pose  interesting  new  problems  in 
structural  dynamics  and  control  In  addition,  modern  sensor  equipment  has  become  more 
and  more  accurate,  making  it  highly  sensitive  to  structural  vibration.  For  example,  an 
antenna  or  optical  sensor,  requiring  a  high  degree  of  pointing  accuracy  may  be  mounted 
on  a  flexible  structure  connected  to  the  spacecraft  bus  The  antenna  will  be  subjected  to 
vibration  from  the  spacecraft  bus  due  to  slew  maneuvers,  momentum  wheel  vibration,  and 
thruster  firing  Passive  vibration  isolation  can  absorb  some  of  the  energy  However, 
modern  payloads  have  led  to  a  requirement  for  active  vibration  control  and  isolation 

There  have  been  two  major  approaches  to  vibration  control  of  flexible  structures. 
The  first,  and  more  conventional  approach,  is  based  on  modal  control  of  the  structure  such 
as  described  by  Newman  [Ref  1]  This  approach  usually  suffers  from  modeling  errors  and 
state  estimation  The  use  of  modal  model  based  control  laws  has  been  an  area  of  intense 
research  during  the  last  decade  The  second  approach  describes  the  structural  response  in 
terms  of  an  elastic  disturbance  which  travels  through  the  structure  [Ref  2]  In  this 
approach,  compensators  are  designed  to  reduce  the  effects  of  incoming  waves  on  outgoing 
waves  The  advantages  of  the  wave  absorbing  approach  are: 

•Relatively  simple  to  implement 

•Broadband  control 

•Does  not  require  finite  element  analysis  or  modal  model 


B.  FOCUS  OF  THESIS 


This  thesis  concentrates  on  three  areas.  First,  a  wave  absorbing  controller  will  be 
derived  using  an  H-infinity  approach  similar  to  one  used  by  Matsuda  and  Fuji  [Ref  3], 
Second,  the  control  law  will  be  simulated  using  a  finite  element  model  of  the  beam  to  be 
studied  The  transfer  function  of  the  open  loop  system  will  be  found  using  modal 
truncation,  then  the  loop  will  be  closed  with  the  H  infinity  controller  Last,  the  control 
law  will  be  implemented  on  the  Naval  Postgraduate  School's  Flexible  Spacecraft  Simulator 
(FSS).  The  control  law  will  be  applied  by  a  real  time  control  system  using  Matrix^™, 
System  Build™,  and  AC- 100  by  Integrated  Systems  Incorporated  (ISI)  The  results  of 
the  experiment  will  be  compared  the  analytical  computer  simulation 
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II.  THEORETICAL  ANALYSIS 


A.  WAVE  ABSORBING  CONTROL 

This  thesis  approaches  active  control  of  the  Flexible  Spacecraft  Simulator  (FSS) 
using  a  wave  absorbing  control  theory.  The  response  of  the  beam  is  viewed  as  a  travelling 
elastic  disturbance  due  to  a  locally  applied  force.  This  approach  has  been  developed  by  a 
number  of  sources  such  as  von  Flotow  and  Shafer  [Ref  4],  MacMartin  and  Hall  [Ref  5], 
Matsuda  and  Fuji  [Ref  3]  and  Fuji  et  al.  [Ref  6], 

1.  Dynamic  Model 

The  FSS  setup  is  shown  in  Figure  1 .  The  base  of  the  flexible  beam  is  fixed  to  the 
center  body  with  the  other  end  of  the  beam  free.  Using  the  Euler-Bernoulli  beam  theory 
the  dynamics  of  a  flexible  beam  are  given  by  the  following  partial  differential  equation: 

=  O  (2.1) 
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The  boundary  conditions  are  : 


Clamped  end: 

w(0,0  =  0,  w(0,0  =  0 

(2.2) 

Actuator  moment: 

EI^0,t)  =  M, 

(2.3) 

where  w(x,t)  is  the  lateral  displacement,  is  the  bending  moment,  El  is  the  bending 
rigidity,  and  pA  is  the  mass  per  unit  length.  The  moment  and  shear  force  at  the  free 
end  are  zero,  but  will  not  be  necessary  for  the  derivation.  As  can  be  seen  in  equation 

(2.3),  the  control  forces  can  be  accounted  for  via  the  boundary  conditions,  which 
allows  the  force  term  on  the  right  hand  side  of  equation  (2.1)  to  be  taken  as  zero, 
a.  Laplace  Transform  Solution  to  Beam  Equation 

The  partial  differential  equation  (2.1)  can  be  analyzed  using  the 

Laplace  transform.  The  initial  conditions,  ^(x,  0)  and  w(x,  0)  are  assumed  to  be  zero. 
The  nonzero  values  of  w(x,  0)  and  w(x,  0)  obviously  do  not  make  a  difference  in  the 

final  result. 


EMx,  0]  =  lV(x,  ■'■)  =  Jr  '•’(x.  t)e-''dl 

(2.4) 

Ll^]  =  x^W(x,.s) 

(2.5) 

(2.6) 

The  fourth  order  partial  derivative  becomes  an  ordinary  derivative  since  time  is 
"frozen"  in  the  s  domain  This  reduces  the  Euler-Bernoulli  beam  equation  to: 


a 


yd^W{x,s)  ,  2 


dx^ 


+  \2^(x,.v)  =  0 


(2.7) 


where 
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(2.8) 


pA 

b.  State  Space  Representation 

The  analysis  of  equation  (2.7)  can  be  simplified  by  using  state 
space  methods  A  state  vector  containing  the  cross  sectional  properties  of  the  beam  is 
introduced. 


where, 


w  =  sw 

dx 

M  =  EI^ 

dx~ 

O  = 


is  the  lateral  velocity 
is  the  angular  velocity 

is  the  internal  bending  moment 
is  the  internal  shear  moment 


The  state  vector  (2  9)  can  be  used  to  transform  equation  (2.7)  to  a  set  of 
first  order  state  equations  in  terms  of  x  in  the  form: 


(2.10) 


Let  us  introduce  p  =  ^  and  proceed  as  follows 


^'■'1  dw  A 

dx  dx  ^ 


(2.11) 


dvi  d-w  A/  .V 

iy=^Tx  =  ^ix^  =  ‘'^¥i  =  -^y^-py^ 


(2.12) 
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(2.13) 


This  leads  to  the  state  equation: 

0  10  0 
0  0/70 
0  0  0  1 
-/7  0  0  0 

The  solution  of  the  ordinary  differential  equation  shown  above  is  given 
by 

r(jc)  =  e^^V(0) 

^Ue'^U-'YiO)  (2.18) 

where  A  is  found  by  similarity  transformation  of  A  .  This  can  be  done  by  solving  the 
following  eigenvalue  problem 


UA  =  AII 


(2.19) 
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(2.20) 


'  X]  0  0  o’ 

A  -  [E  0  Xz  0  0 

■>  2  0  0  ?L3  0 

0  0  0  ^4 

’  1  1  1  1  ’ 

X,|  Xz  X3  X4 

u=  ^  ^  ^  ^  (2.21) 

p  p  p  p 

X.J  xl  >.4 

_T  T  T  T  _ 

where  the  eigenvalues  are 

1+/  ' 

X=\  i  (2.22) 

-1  +  / 

-1-/ 

As  can  be  seen  from  equations  (2. 1 8-22),  the  solution  yields  imaginary 
solutions.  A  more  useful  solution  can  be  transforming  the  matrix  A  or  its  diagonal 
transformation  to  real  Jordan  canonical  form.  [Ref  3]  As  discussed  in  the  next  section 
this  will  also  yield  a  state  vector  in  terms  of  travelling  waves  instead  of  cross  sectional 
variables 

c.  State  Vector  in  Terms  of  Travelling  Waves 

The  state  vector  Y  from  the  previous  section  is  defined  in  terms  of 
the  cross  sectional  variables  of  the  beam.  As  stated  at  the  beginning  of  this  chapter,  we 
desire  to  study  the  disturbances  in  the  beam  in  terms  of  travelling  waves.  The  matrix  A 
can  be  transformed  to  real  Jordan  canonical  form  by  the  matrix  K  [Ref  7,  pp  144,145] 
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i-jo  0 

i  i  0  0 

0  0  1-5 

0  0  i  5 


(2.23) 


The  new  matrix  A  can  then  be  formed  by  the  similarity  transformation 


A  =  iK-W-^)AiUK)  =  P-^AP 


(2.24) 


where 


VT/7  0  Jlp  0 


P=IJK  =  ^ 


M  — r)3 


0  Jlp  0  -i/2p 

T  T  -)  2 

-/?" 


(2.25) 


Using  this  transformati  :  n  the  state  vector  can  be  changed  to 


Y  =  PV 


(2.26) 


where  the  new  state  vector  V  contains  the  amplitudes  of  the  travelling  wave  modes. 


(2.27) 


This  notation  is  similar  to  that  used  in  liturature  used  in  microwave  circuits.  [Ref  4,  pg. 
674]  At  any  point  along  the  beam  the  motion  can  be  described  by  four  modes,  two 
incoming  and  two  outgoing  waves  as  shown  in  Figure  2. 
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Figure  2.  Travelling  Wave  Modes 


The  new  state  vector  varies  with  postion  along  the  beam  according  to  the 


new  state  equation 


110  0 

110  0 

0  0-1  1 

0  0-1-1 


(2.28) 


Equation  (2  28)  has  a  solution  of  the  form 


V(x)  =  e^'^V(0) 


(2.29) 


The  matrix  can  be  represented  by 


Q 

0  6-'-' 


(2.30) 


where  A  \  and  A  2  are  the  block  diagonals  of  A 


Next  the  inverse  Laplace  transform  (with  respect  to  x) 

e^=L-^[{sI-Ay^]  (2.31) 

is  used  to  find  the  exponentials  along  the  diagonal  of  equation  (2.30).  Let  C  =  , 

then 


Cv/-/I,)  = 


(2.32) 


(.V/ 


(2.33) 


By  taking  the  inverse  Laplace  transform  using  equation  (2  31)  the  first  diagonal 
element  of  equation  (2.30)  is  found 

,  e^^coscx 

^Ai-x  _ 

-e^sincjc  e^^coscx 

The  second  diagonal  element  can  be  found  in  a  similar  fashion  leading  to  the  solution 

e^'^coscx  e^sincx  0  0 

£;"cos6X  0  0 

0  0  <?“"cos6X  e'‘‘^*sincr 

0  0  e~^^coscx 

To  determine  the  wave  amplitudes  for  the  pinned  and  free  end  of  the  beam  in  Figure 
2,  the  solution  is  determined  for  x  L  This  leads  to  the  same  solution  as  determined 
by  Matsuda  and  Fuji  [Ref  3] 


(2.34) 
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e^coscjc  e'^^smcx 
-e^^sincx  e^coscx 


e  ‘^■^coscjc  e  “^^sincjc 
-e~‘'^smcx  £?”“coscx 


(2.36) 


d  Boundary  Conditions 

The  boundary  conditions  are  given  by  equations  (2.2)  and  (2.3).  The 
boundary  conditions  can  be  transformed  to  the  new  coordinates  by  equation  (2.25) 


10  0  0 

0  0  f  0 


(2.37) 


which  can  be  rewritten  as 


1  0  0  0  1^, 

0  0  f  0 


Jlp  0  -Jlp 

p3  -p3 

0  0 


-p3  p3 


-Jlp 


(2.38) 


{Jlp  0  {Jip  0 
0  %j2p  0  -fMp 


(2.39) 


This  equation  can  be  solved  for  the  outgoing  wave  amplitudes  [h, ,  />,]  as  a  function  of 
the  incoming  wave  modes  [  ,«,]  and  the  control  moment 
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-1  0 
0  1 


a\ 

ai 


+ 


(2.40) 


2.  Control  Law 

a.  System  Equation  in  Four  Block  Form 

The  wave  equations  derived  in  Section  1  can  be  written  as  a 
standard  four  block  form  in  robust  control  theory.  The  equations  are  as  follows. 


h  =  Sa->r  Bu 

(2.41) 

y  =Ta  +  Gu 

(2.42) 

where  the  matrix  .V,  which  in  this  special  case  is  called  the  scattering  matrix  contains 
the  reflection  coefficients.  The  standard  four  block  form  is  written  as 


(2.43) 


The  block  diagram  for  this  system  is  shown  in  Figure  3 


Figure  3  System  Block  Diagram 
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The  controller  inputs  are  of  the  form 


u  =  Ky  (2.44) 

tMc  =  KQ 

El 


Combining  compensator  K  with  the  output  equation  given  by  (2.42)  gives 


u  =  KTa  +  KGu 

(2.45) 

or 

u  =  HTa 

(2.46) 

where 

II 

1 

1 

(2.47) 

or,  in  other  words. 

K  =  H{I  +  HCjy^ 

(2.48) 

The  closed  loop  system  between  a  and  h  can  be  determined  by  combining  equation  (2.46) 
with  (2.41) 


h  =  Sa  +  BHTa  =  {S  +  BHT)a  (2.49) 

or 

h  =  Scia  (2.50) 


where  .S'^,  is  the  closed  loop  scattering  matrix 

For  the  system  derived  in  Section  1 ,  the  four  block  form  equation  can  be 
determined  from  the  transformation  (2.26)  and  equation  (2.40). 


w 

Jlp  0  Jlp  0 

^  1  1  *7 

a\ 

e 

s 

p^^  pt  -pi  pi 

\2j 

1 

o 

o 

b^ 

1  f  y  1 

"in")  n 

-p'^  p~^  p~^  p~^ 

, 
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The  input  equation  (2.44)  can  be  determined  by  using  the  second  and  third  row  equations 
in  (2.51)  above. 

Q  =  \p^a\+\p^a2-\ph\+\ph2  (2.52) 

fjMc  =  J2pa2  -  J2pb2  (2.53) 

The  above  two  equations  are  combined  producing 


0  =  p^-a\  +p-a2  - 


(2.54) 


Equations  (2  40)  and  (2  54)  form  an  equation  appropriate  for  a  control  system  design. 

-10  0 
0  1  -yflp-^ 

-(f)  =  _ 

The  first  row  of  this  equation  can  be  dropped  since  cannot  be  physically  controlled  by 
M  This  wave  mode  is  called  a  "far  field"  mode  since  it  represents  a  mode  transmitted 

c 

along  the  beam's  axis  for  a  long  distance.  Therefore,  the  closed  loop  system  S^,  can  be 
determined. 

.S'rf  =  [o  1  ]  +  (-V2p-')/{p=  p\  ] 

^[-X\-X]  (2-56) 

where 

JlHp'^-  (2.57) 


a\ 

<  a2  "  (2.55) 

Me  \ 
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b.  H  infinity  Compensator  Design 

In  order  to  minimize  the  effect  of  the  incoming  waves  on  the  outgoing 
waves,  it  is  desirable  to  choose  a  compensator  K  for  which  the  "magnitude"  of  the  closed 
loop  scattering  matrix  is  minimized  This  is  motivated  by  drawing  energy  out  of  the 
system  at  the  junction  point  where  the  actuator  is  located. 

The  magnitude  of  a  matrix  can  be  measured  by  its  norm.  One  type  of 
matrix  norm  is  the  H  infinity  norm  which  is  given  by 

\\Sci\L  =  Sup(G[Sci(J'^)])  (2.58) 

G  is  the  largest  singular  value  The  singular  values  can  be  found  by 

cfSci)  =  jMSliSci)  (2.59) 

where  is  the  complex  conjugate  of  Sch  and  A,,  are  the  eigenvalues  of  the  matrix  in 
parentheses.  The  eigenvalues  are  found  by 


-X(-jrn) 

1  -X(-jrn) 


-X(jr^)  I  -X(;ra)  ] 


(2.60) 


A  -  X(jw)X{-jvi)  X{-im)[  1  -  Xijvs)] 
A^(/m)[l -X{-im)]  A-  [1  -X(jm)][\  -X(-/tn)] 


Solving  the  determinant  above  leads  to 

=  X  =  2[x(/ra)  -  j][ji'(-/ro)  -  j]  +  j  (2.«2) 

The  smallest  solution  for  equation  (2.62)  occurs  at  A1(.S')  =  j  which  gives  G  =  Now 
equation  (2  57)  can  be  solved  for  H. 
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(2.63) 


This  transfer  function  was  derived  for  a  beam  clamped  at  the  base  and  free 
at  the  end  The  transfer  function  above  was  also  derived  for  a  pinned-free  beam  by 
Matsuda  and  Fuji  [Ref  3]  using  the  H-infmity  approach  A  similar  tranfer  function  was 
found  by  MacMartin  and  Hall  [Ref  5]  by  minimizing  the  H  infinity  norm  of  the  power 
flow  into  the  system  MacMartin  and  Hah  uiso  point  out  that  this  solution  extracts  half  of 
the  power  at  all  frequencies  This  is  in  contrast  to  other  controllers  which  are  more 
effective  at  a  given  mode,  but  cannot  add  significant  damping  at  all  of  the  modes  at  the 
same  time 

The  control  law  given  in  (2.66)  is  effectively  a  "half  differentiator"  It  is 

similar  to  velocity  feedback  with  a  forty  five  instead  of  a  ninety  degree  phase  lead  In 
order  to  implement  this  control  law  the  function  Js  must  be  estimated  by  a  rational 
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transfer  function.  As  shown  by  MacMartin  and  Hall  [Ref  5],  this  irrational  function  can 
be  approximated  over  a  wide  frequency  range  by  a  finite  number  of  zeros  and  poles  spaced 
logrithmically  along  the  negative  real  axis.  Figure  4  shows  the  Bode  plot  for  four  zeros 
and  four  poles. 


(i+io-^)  (j+io--)  (^+10°)  (j+io-) 

(j+io-^)  (j+10-')  (j+10')  (j+10^) 


(2.67) 


Figure  4.  Rational  Estimate  of  Js 
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B.  SENSOR  AND  ACTUATOR 

The  sensors  and  actuators  for  the  FSS  beam  are  nearly  collocated  at  the  base  and 
can  be  seen  in  Figure  5.  Both  the  sensor  and  actuator  are  made  of  piezoelectric  ceramic 
As  seen  in  the  figure,  the  piezoelectric  wafers  are  stacked  This  allows  for  high  sensitivity 
from  the  sensor  and  more  actuator  power  in  a  small  area. 


Figure  5.  First  Beam  Element  with  Sensor  Actuator  Pair 


1.  Piezoelectric  Material 

Piezoelectric  ceramics  have  many  advantages  as  actuators  and  sensors.  Table  1 
shows  the  advantages  of  using  piezoelectric  materials  in  structural  vibration  control. 

[  Ref  8] 

The  usefiilness  of  piezoelectric  ceramics  in  vibration  control  is  derived  from  their 
ability  to  convert  electrical  energy  into  mechanical  energy  and  vice  versa  When  a  force  is 
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applied  to  the  material,  it  creates  a  voltage  proportional  to  the  applied  force  Conversly, 
when  an  electric  potential  is  applied  across  the  ceramic,  it  creates  mechanical  strain. 


Sensor 

Actuator 

high  strain  sensitivity 

low  noise 

low-moderate  tempurature  sensitivity 
easily  implemented 

high  stiffness 

sufficient  stress  to  control  vibrations 

good  linearity 

low  power  consumption 

easily  implemented 

Table  1.  Advantages  of  Piezoelectric  Ceramics 


2.  Sensor 

It  can  be  shown  that  a  piezoelectric  material  produces  voltage  proportional  to 
strain  when  force  is  applied  [Ref  9]  This  voltage  is  given  by  the  equation 

Kv  =  /.v(^)(l -v)e.,  (2.68) 


The  material  constants  for  the  Navy  Type  II  PZT,  which  is  used  in  this  study,  are  given  in 
Table  2. 


Constant 

Description 

Units 

Value 

d3, 

Lateral  Strain 

m/V  or  Coul/N 

1.8e-10 

E 

Young's  Modulus 

Pascal  (N/m“) 

6.30el0 

V 

Poissons  Ratio 

n/a 

0.35 

D 

Abs  Permittivity 

Farad/m  or  N/V' 

1.5e-8 

ts 

Sensor  Thickness 

m 

2  X  1 .905e-4 

tb 

Beam  Thickness 

m 

1  588e-3 

A 

Sensor  Area 

m- 

4.839e-4 

Table  2  Material  Constants  for  Navy  Type  II  PZT 
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The  strain  developed  at  the  cross  section  of  the  beam  at  a  distance  h  from  the  midplane  of 
the  beam  is 


(2.69) 


Combine  this  equation  with  (2.68)  to  get 

V,(x)  =  (2.70) 

where 

k  =  -ht,^{\-v)  (2.71) 

Since  the  sensor  is  distributed  along  the  beam,  however,  the  measured  sensor  output  is 
the  average  voltage  over  the  length,  L,  of  the  wafer. 


(2.72) 


(2.73) 


Thus,  the  voltage  output  of  the  sensor  is  directly  proportional  to  the  effective  bending 
angle  over  the  sensor  element.  The  value  of  ^  for  the  FSS  beam  is  -1.378  x  10"^ 
rad/volt 

3.  Actuator 

The  actuator  produces  a  moment  proportional  to  the  voltage  across  its  terminals 
The  electric  field  across  the  wafer  is  related  to  the  voltage  by 
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(2.74) 


Thus,  the  actuator  produces  a  moment  which  is  proportional  to  the  voltage  applied  across 
its  terminals 
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III.  STRUCTURAL  RESPONSE  SIMULATION 


The  open  loop  transfer  function  was  determined  by  a  modal  model  which  was 
generated  by  truncating  a  desired  number  of  modes  The  open  loop  transfer  function,  G(s), 
and  the  H-infmity  controller,  //(.?),  were  then  discretized  and  formed  into  the  closed  loop 
system  as  shown  in  Figure  6 


reference 
input 


Beam 

Modal  Model 
G(z) 


theta(z) 

- ^ - ► 


Controller 

H(z) 


*4 


Figure  6.  System  Block  Diagram 


A.  OPEN  LOOP  TRANSFER  FUNCTION 

The  beam  dynamics  were  determined  by  a  modal  model  The  equation  of  motion 
for  a  lumped  parameter  model  is  given  by  the  matrix  equation 


Mq  +  Cq  +  Kq  =  F(t) 


(3.1) 


where 

M  is  the  mass  matrix 
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C  is  the  damping  matrix 
K  is  the  stiffness  matrix 
q  is  the  vector  of  generalized  coordinates 

The  beam  is  modeled  with  eight  finite  elements  as  shown  in  Figure  7  The  mass 
and  stiffness  matrices  were  constructed  using  the  finite  element  method.  [Ref  10,  pp. 
300-328]  Element  mass  and  stiffness  matrices  were  calculated  and  combined  to  form  the 
global  mass  and  stiffness  matrices  for  equation  (3.1).  These  matrices  can  be  seen  in  the 
program  BEAM.M  in  the  Appendix  A  It  should  be  noted  that  the  vertical  dispacements 
and  rotation  at  the  base  are  equal  to  zero  due  to  clamped  boundary  conditions. 
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The  generalized  coordinate  vector  can  be  decoupled  into  modal  coordinates.  [Ref.  11, 
pp.  104- 107] 

Wi 

01 

^0)  =  ]  *  i  =  OriCO  =  t  ({),  11,(0  (3.2) 

•  /=i 

-  08  , 

The  vector  r|(0  contains  the  modal  coordinates,  and  the  modal  matrix  O  contains  the 
system  eigenvectors  normalized  with  respect  to  the  mass  matrix.  The  number  of  modes, 

«,  for  the  lumped  parameter  model  is  equal  to  the  number  of  generalized  coordinates.  A 
reasonable  solution  can  be  found  by  truncating  the  number  of  modes  into  a  physically 
meaningful  number  For  example,  three  modes  are  often  adequate  to  describe  the 
longitudinal  vibration  of  a  simple  cantilevered  beam.  [Ref  11,  pp.  292-293] 

Equation  (3.2)  is  substituted  for  the  generalized  coordinates  in  equation  (3  1)  and 
premultiplied  by  to  get  the  equation  of  motion  in  terms  of  the  modal  coordinates 

fi-i-A,  ii-i- A,tr|  =  (3.3) 

The  first  coefficient  matrix  is  the  modal  damping  matrix 

IpiTHi  0  0  0  0 

0  .00  0 

A,.  =  cD^(’cD=  0  0  *  0  0  (3.4) 

0  0  0  *  0 

0  0  0  0  2pi6TiT|6 
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The  modal  stiffness  matrix  is 


~  -nj?  0  0  0  0  ' 

0  TiJ2  0  0  0 

=  0  0  *  0  0  (3.5) 

0  0  0  *0 

0  0  0  0 

and  the  modal  force  vector  is 

f  '' 

0 

Me 

0 

F(/)  =  O^F(0  =  OM  *  \  (3.6) 

• 

0 

0 

Note  that  since  the  actuator  is  mounted  on  the  first  element,  the  external  moment.  A/  ,  is 

applied  to  the  second  coordinate,  q{2)  =  0 1 

The  modal  matix  equation  (3.3)  can  be  separated  into  its  individual  equations  in 

the  form 

f),  +  2p,Ti5,f)/  +  tu^ri,  =fiU{t)  i  I,2,...,n  (3.7) 

The  Laplace  transform  of  this  equation  is. 

s^r\ ,  +  .s^p/tn,  ri  /(.v)  +  ,(.s)  =  / ,  (3-8) 

which  leads  to  the  transfer  function  for  q,(.s)  and  (I(s)  as  follows 
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ri,(-y)  _  fi 

s^+lpWjS+xnj 


(3.9) 


Equation  (3.9)  gives  the  transfer  function  from  the  input,  U(s),  to  the  i-th  modal 
coordinate,!)/,  where  /  =  1,2,  .  ,16. 

Now,  the  system  transfer  function  for  a  physical  coordinate  can  be  determined. 
The  output  equation  for  a  single  input,  multiple  output  system  is 

yit)  =  C,q{t)  (3.10) 

As  discussed  in  Chapter  II,  the  piezoelectric  sensor  mounted  on  the  first  element  measures 
the  angle  0 1  Thus  the  output  equation  becomes 

y(t)  =  0,(0  =  [  0  1  0  •  •  0  ]q{t)  (3.11) 

0i(O  =  ^2(O  (3.12) 

This  equation  leads  to  a  single  input,  single  output  system.  Using  the  transformation  to 
modal  coordinates  given  by  equation  (3.2),  the  output  equation  becomes 

0,(0  =i  (1)2,, Tl/(0  (3.13) 

1=1 

where  <1)2.,  denotes  elements  of  second  row  of  O  Convert  this  equation  to  the  s  domain 
using  the  Laplace  transform  Therefore, 

e  I  (.v)  =  S,|>2./n  .(■'■)  (s-'-t) 

,=  l 

Combining  this  equation  with  (3  9),  and  the  fact  that  the  input  lJ(s)  M/s),  yields  the  open 
loop  transfer  function 


(3.15) 


ei(^)  _ 

Mc{s)  y_j  s^+lpimjS+xn^ 

As  seen  in  Figure  8,  one  way  to  view  equation  (3.15)  is  n  single  transfer  functions 
combined  in  parallel.  This  can  be  easily  implemented  using  the  Matlab™  Control  Systems 
Toolbox  "parallel. m"  command 


Figure  8  Open  Loop  Transfer  Function 


B.  DISCRETE  TRANSFER  FUNCTIONS 

Each  transfer  function  was  tranformed  to  the  z  (discrete)  domain  before  forming 
the  closed  loop  system  The  relationship  between  the  s  (continuous)  and  z  (discrete),  is 

2  =  (3.16) 

where  T  is  the  sampling  period  and  s  =  —a  +  jxnl 

The  open  loop  transfer  function,  determined  by  equation  (3  15)  was  converted  to  a 

discrete  transfer  function  using  a  Tustin  transformation  The  Tustin  transformation  shown 

below  is  the  Fade  approximation  to  the  exponential  (3.16).  [Ref  12,  pp.  253-282] 
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(3.17) 


The  Tustin  transformation  maintains  the  frequency  response  of  the  continuous  system 
while  preserving  the  mapping  of  the  s-plane  into  the  z-plane. 

The  transfer  function  for  the  estimated  H-infinity  controller  was  transformed  to  the 
z  domain  by  the  matched  pole  zero  technique.  Recall  from  Chapter  II  that  the  H-infmity 
controller  was  irrational  and  was  estimated  using  zeros  and  poles  spaced  logrithmically 
along  the  negative  real  axis.  These  zeros  and  poles  can  be  mapped  directly  to  the  z 

domain  by  replacing  each  zero  or  pole  term  5  +  a  with  its  discrete  equivalent,  z  - 
The  gain  for  the  discrete  transfer  function  is  chosen  so  that  G{z)\z=\  =  G(5)l  j=o . 

[Ref  12,  pp.  3  04-3 05]  Using  Matlab™  the  equivalent  form  of  equation  (2.64)  is 


^ _  ,(2+ 1  ■  000000)(z+. 99995 )(z+. 995 0 1 2)(z+. 6065 3 1 ) 
-  ^(z+.999995)(z+.999500)(z+.951229)(z+.006738) 


(3.18) 


C.  CLOSED  LOOP  SYSTEM 

Control  of  the  beam  was  simulated  by  closing  the  loop  using  H  infinity  controller, 
(yj )  and  derivative  (.v).  Also,  a  second  order  filter  was  incorporated  to  study  the  effects 
of  filtering  the  sensor  output  in  the  system  The  closed  loop  transfer  function  was 

determined  using  the  feedback  command  in  Matlab''^. 

D.  SIMULATION  RESULTS 

The  open  and  closed  loop  Bode  plots  for  the  H-infinity  controller  are  shown  in 
Figure  9.  The  plot  shows  a  reduction  in  magnitude  over  a  broad  frequency  range  as 
expected. 
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Figure  9  Open  and  Closed  Loop  Bode  Plots 
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The  Bode  plot  becomes  inaccurate  at  frequencies  above  the  Nyquist  rate  of  3 14  radians 
per  second  due  to  the  system  being  discrete  with  a  sampling  period  of  0.01  seconds. 

Figure  10  shows  the  time  response  to  a  sinusoidal  input  of  1  Hz  which  is  the  first 
natural  frequency  of  the  beam  This  plot  shows  a  dramatic  reduction  in  the  amplitude  of 
the  closed  loop  system  for  the  first  two  modes  ( 1  and  1 7  Hz) 
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IV.  EXPERIMENTAL  ANALYSIS 


A.  SET  UP 

The  wave  absorbing  controller  was  implemented  on  the  Flexible  Spacecraft 
Simulator  at  the  Naval  Postgraduate  School  Figure  1 1  shows  an  overall  picture  of  the 
system. 


Figure  1 1 .  Flexible  Spacecraft  Simulator 


The  spacecraft  simulator  consists  of  a  center  body  which  floats  over  a  granite  table 
on  air  pads  and  rotates  around  an  air  bearing  A  momentum  wheel  is  mounted  on  the  body 
to  allow  for  slew  maneuvers  and  in  this  experiment  is  used  as  a  disturbance  source  The 
flexible  beam  is  attached  to  the  center  body  as  shown  in  Figure  12. 


Figure  12.  Flexible  Beam 

The  beam  is  controlled  by  a  piezoelectric  sensor-actuator  pair  mounted  at  the  base  of  the 
beam  as  seen  in  Figure  13.  The  dimensions  of  the  sensor  and  actuator  are  shown  in  Figure 
5,  Chapter  II. 


Figure  1 3  Sensor  and  Actuator 


The  digital  control  system  was  implemented  using  Matrix/*^,  System  Build™,  and 
AC- 100™  control  hardware.  The  control  law  was  designed  using  block  diagrams  in 
System  Build™'  The  system  was  then  compiled,  linked  and  loaded  into  the  AC- 100™,  a 
real  time  controller  Figure  14  shows  the  control  software  block  diagram 
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Discrete  Super-Block  Sampling  Interval  First  San^le  Ext. Inputs  Ext. Outputs  Enable 

visl  0.0100  0,  18  9  Parent 


Figure  14  Control  System  Block  Diagram 
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O.OOSs 


The  AC- 100™  has  an  output  of  plus  or  minus  10  Volts.  This  signal  was  amplified 
by  two  Kepco™  amplifiers  connected  in  series.  This  allowed  a  range  of  plus  or  minus  75 
Volts  for  the  piezoelectric  actuator  input. 

The  momentum  wheel  was  used  as  a  disturbance  force  and  is  controlled  by  voltage 
with  1  Volt  per  300  RPM.  In  order  to  provide  a  disturbance  torque  at  a  desired  frequency 
the  momentum  wheel  was  spun  up  to  a  set  RPM,  then  the  wheel  RPM  was  modulated  at 
the  desired  frequency.  The  control  blocks  for  the  momentum  wheel  are  also  shown  in 
Figure  14. 

B.  EXPERIMENTAL  RESULTS 

1.  Initial  Tests 

The  initial  choice  for  the  digital  sampling  rate  was  100  Hz.  Sampling  rates  of  50, 
200,  and  400  Hz.  were  also  used.  The  100  and  200  Hz  sampling  rates  appeared  to  be  a 
reasonable  choice. 

The  main  difficulty  at  the  initial  stages  of  the  experiment  was  high  frequency  noise 
coming  from  the  sensor.  Whenever  there  was  a  large  jump  in  the  sensor  signal  between 
samples,  the  output  to  the  actuator  would  have  a  large  jump  in  voltage,  causing  the 
actuator  to  click.  These  voltage  jumps,  possibly  as  high  as  150  Volts,  could  cause  damage 
to  the  piezoelectic  actuator.  Thus,  different  sampling  rates  were  tried.  The  higher 
sampling  rates  worked  well  at  reducing  the  voltage  steps  between  samples,  since  the 
sampling  period  was  shorter,  making  the  actuator  output  closer  to  continous.  However, 
the  higher  sampling  rates  also  allowed  more  of  the  higher  frequency  noise  into  the  system 
The  lower  sampling  rates,  allowed  less  noise  in  the  system,  but  had  large  voltage  jumps 
between  samples. 

The  solution  to  these  problems  was  to  insert  a  low  pass  filter  at  the  sensor  output. 
A  number  of  filters  of  different  orders  and  cutoff  frequencies  were  tried  The  final  choice 
was  a  second  order  filter  with  a  damping  ratio  of  0.707  and  a  cutoff  frequency  of  30 
radians  per  second  (4.8  Hz).  This  filter,  combined  with  a  100  Hz  sampling  rate,  allowed 
good,  stable  performance  for  H-infmity  and  derivative  control. 
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2.  Free  Response  of  the  Beam 

The  free  response  of  the  beam  was  tested  for  the  first  two  modes.  The  beam  was 
set  to  an  initial  condition  and  allowed  to  vibrate  freely  The  plots  for  response  of  the  first 
two  modes  are  shown  in  Figures  15  and  16.  Note  that  the  second  mode  plot  data  was 
collected  at  a  rate  of  500  Hz  for  more  accuracy.  In  order  to  get  a  good  plot  of  the  second 
mode  the  data  was  filtered  using  a  high  pass  FIR  filter  with  a  cutoff  frequency  of  14  Hz 
This  filtered  out  the  dominant  first  mode  and  allowed  for  a  better  view  of  the  second 
mode. 

The  damping  ratio  of  each  mode  was  determined  by  noting  that  the  amplitude  of 
the  wave  at  any  time,  t,  is  determined  by 

A(t)  =  (4.1) 


where  A„  is  the  initial  amplitude  and  is  the  frequency  of  that  mode  If  the  amplitude  is 
measured  at  two  different  points  on  the  plot  the  damping  ratio  can  be  found  by  solving  for 


(4.2) 


The  damping  ratios  of  the  first  two  modes  are  shown  in  Table  3. 


Mode 

Frequency  (Hz  /  (rad/sec) ) 

Damping  Ratio 

1 

0.95  /  5.97 

3.7x10-' 

2 

17.0/  106.8 

3.9x10' 

Table  3  Free  response  damping  ratios 


The  damping  ratios  determined  experimentally  were  used  for  the  finite  element 
program  described  in  Chapter  III  The  higher  modes  were  estimated  to  be  about  the  same 
for  the  program.  Also,  note  that  the  first  two  natural  frequencies  agree  well  with  the 
program  frequencies  of  1.01  and  17  8  Hz 
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3.  Initial  Condition  Response  with  Active  Control 

In  this  part  of  the  experiment,  the  response  of  the  beam  to  an  initial  condition  was 

tested  with  two  different  control  laws.  The  first  was  the  H-infinity  (or  ) 
controller,  and  the  second  was  a  derivative  control.  As  pointed  out  in  Chapter  II,  the 
H-inifity  controller  in  this  case  is  similar  to  the  derivative  control.  It  provides  a  phase  lead 
of  45  degrees  instead  of  90  for  the  derivative  controller  Its  magnitude  rolls  off  more 
sharply  compared  to  the  derivative  controller  at  higher  frequencies. 

The  results  are  shown  in  Figures  17  and  18.  A  damping  ratio  for  each  case  was 
calculated  in  a  similar  fashion  to  the  free  response  in  the  previous  section.  The  measured 
frequencies  and  calculated  damping  ratios  are  shown  in  Table  4. 


Frequency  Damping 

(Hz) _ Ratio 


Free  Response  0.95  3.9x10'^ 

H-infinity  Control  0.98  3.7x10  ' 

Derivative  Control  0.96  3.9x10" 


Table  4  Damping  Ratios 

From  the  results  in  Table  4  the  derivative  controller  seems  to  be  slightly  better  than 
the  H-infinity  controller  at  this  frequency.  However,  as  can  be  seen  in  Figure  18,  the 
derivative  controller  began  to  excite  the  second  mode  after  about  30  seconds.  The 
derivative  controller  was  also  much  more  sensitive  to  noise  in  the  system  and  often  caused 
spikes  in  the  voltage  signal  sent  to  the  actuator.  The  H-infinty  controller,  on  the  other 
hand,  was  less  sensitive  to  noise  and  gave  a  smoother  input  to  the  actuator  The  H-infinity 
controller  turns  out  to  produce  reliable  performance  and  was  more  stable 


N (stcxlOO) 

B«am  Angle  ~  H^nftnlty  Control  {T«.01) 

Figure  17.  Initial  Condition  Response  using  H-infinity  Control 


Figure  18  Initial  Condition  Response  using  Derivative  Control 


4.  Forced  Response 

The  forced  response  of  the  closed  loop  system  was  tested  using  the  momentum 
wheel  as  a  disturbance  source  Different  runs  were  made  using  different  frequencies. 
Both  H-infmity  and  derivative  control  laws  were  implemented 

As  discussed  previously,  the  system  had  stability  problems  with  high  frequency 
noise  and  large  voltage  jumps  Thus,  there  was  a  necessity  to  use  a  low  pass  filter. 
Initially  first  and  second  order  filters  with  30,  50  and  100  rad/sec  cutoff  frequencies  were 
used.  The  most  stable  performance  was  achieved  using  a  second  order  30  rad/sec  filter. 
This  reduced  the  effective  bandwidth  of  the  controller  to  4.8  Hz  Also,  the  filter  began  to 
reduce  the  amount  of  phase  lead  of  the  controller  as  the  disturbance  frequency  came  near 
4  8  Hz 

The  controller  did,  however,  work  well  within  the  low  pass  bandwidth  with 
significant  active  control  action  near  the  first  natural  mode  of  1  Hz.  The  results  of 
closing  the  control  loop  are  given  in  Figures  19  and  20.  The  decrease  in  amplitude  for 
both  controllers  was  calculated  and  is  shown  in  Table  5 


Controller 

(Initial/Final)  amplitude 

Decrease  in  amplitude(dB) 

H-infinity 

5 

14 

Derivative 

3.35 

10.5 

Table  5  Controller  performance  (1  Hz  disturbance) 


The  H-infinity  controller  performed  better  in  this  case  A  significant  reason  for 
this  was  that  the  H-infinty  controller  was  less  sensitive  to  noise  and,  therefore,  the  gain 
could  be  increased  to  get  full  use  of  the  150  V  peak  to  peak  output  of  the  actuator 
amplifier  The  derivative  controller  became  unstable  if  the  gain  was  increased  too  much, 
even  with  the  low  pass  filter  Figure  20  shows  that  the  derivative  controller  begins  to 
excite  higher  modes  after  about  55  seconds  The  H-infinity  controller,  combined  with  the 
second  order  filter,  showed  no  tendency  to  excite  the  higher  modes  when  subjected  to  the 
disturbance  The  H-infmity  controller  had  the  best  performance  with  a  forced  response 
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V.  SUMMARY  AND  CONCLUSIONS 


Active  control  of  flexible  spacecraft  structures  has  become  an  important  topic. 
Larger  and  lighter  structures  supporting  more  sensitive  equipment  have  led  to  the  need  for 
active  control  of  spacecraft  structures  This  thesis  studied  a  wave  absorbing  approach  to 
this  control  problem  A  closed  loop  scattering  matrix  was  derived  which  determines  the 
relationship  between  incoming  waves,  outgoing  waves,  sensor  and  actuator.  The  control 
law  was  found  by  minimizing  the  magnitude  of  this  closed  loop  scattering  matrix. 
Theoretically,  this  approach  allows  for  a  relatively  simple  control  law,  broadband  control, 
and  does  not  require  a  modal  model  of  the  system  for  operation. 

This  control  law  was  then  applied  to  the  Flexible  Spacecraft  Simulator  (FSS).  The 
sensing  and  actuation  for  the  flexible  beam  was  implemented  by  piezoelectric  ceramic 
wafers  mounted  at  the  base  of  the  beam  The  control  law  was  tested  by  simulation  on  a 
computer,  and  then  by  applying  it  to  the  hardware  in  the  lab. 

The  computer  simulation  showed  that  the  controller  was  useful  over  a  broad  range 
of  frequencies  below  the  Nyquist  rate  of  the  digital  system.  The  results  from  the  FSS 
confirmed  this  for  frequencies  below  5  Hz.  At  higher  frequencies  the  controller  had  some 
problems  due  to  sensor  noise  requiring  a  low  pass  filter,  and  thus  limiting  the  bandwidth  of 
the  system  At  the  first  natural  frequency  of  1  Hz  the  H-infinity  controller  performed  well 
and  increased  the  damping  ratio  of  the  system  from  3.93  x  10^  to  3.7  x  10  ^  The 
controller  also  decreased  the  amplitude  due  to  a  1  Hz  sinusoidal  disturbance  by  14  dB. 

Some  recommendations  for  further  study  include  adding  a  weighting  function  to 
the  terms  in  the  closed  loop  scattering  matrix.  This  would  allow  for  weighting  the  output 
of  the  controller  at  desired  frequencies.  It  also  allows  for  implementation  of  a  low  pass 
filter  before  the  H-infinity  norm  of  the  matrix  is  minimized  Another  area  of  study  would 
be  the  use  of  analog  filters  on  either  the  sensor  output  or  the  actuator  input.  The  former 
could  be  used  to  filter  out  noise  from  the  sensor  The  latter  would  be  useful  in  smoothing 
out  the  signal  sent  from  the  power  amplifiers  to  the  piezoelectric  actuators.  The  filters. 
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however,  will  change  the  phase  of  the  closed  loop  system,  which  is  another  area  to  be 
addressed.  Lastly,  the  sensor  and  actuator  were  assumed  to  be  collocated  for  this  study. 
Their  centers  are  1 .96  inches  apart  This  should  not  be  a  problem  for  the  lower  modes  of 
the  system  However,  it  is  possible  that  at  higher  frequencies,  the  phase  difference 
between  the  wave  at  the  sensor  and  actuator  could  cause  a  stability  problem.  This  could 
be  alleviated  by  solving  the  scattering  matrix  for  a  noncollocated  sensor  and  actuator 

In  conclusion,  the  use  of  a  wave  absorbing  controller  with  piezoelectric  sensors 
and  actuators  could  be  a  viable  approach  to  vibration  suppression  of  space  structures  As 
discussed  by  von  Flowtow  [Ref  13]  and  Fuji,  Ohtsuka  and  Murayama  [Ref  6], 
eventually,  this  approach  can  be  applied  to  large  flexible  spacecraft  structures  with 
multiple  noncollocated  sensors  and  actuators.  The  wave  absorbing  control  approach 
combined  with  new  sensor  and  actuator  technology  such  as  piezoelectric  ceramics 
presents  a  wide  area  for  further  research 
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!  Beam  Simulator 


I  This  program  models  a  beam  27  inches  long  with  lumped  masses 
%  to  change  the  natural  vibration  frequencies.  The  beam  is  ^ 

%  modeled  using  a  finite  element  model  with  8  elements.  The 
^  transfer  function  (theta/moment  at  element  1)  of  the  beam  is 
%  found  using  nodal  truncation.  Tha  program  f  ® . 

%  control  laws  to  be  tested  by  the  user.  Included  are  the  H  ^ 

%  infinity  controller  {sqrt(s)),  Derivative,  and  Integral  ^ 

%  control.  Also  a  second  order  filter  with  a  variable  cutoff  ^ 

%  frequency  can  be  applied  to  the  H  infinity  or  Derivative  ^ 

%  controller.  * 

J***************************************************************** 

clear 


m=. 112014; 

ml=.4672; 

Il=2.0436e-4; 

h=.0635; 

hl=.1143; 

h8=. 14478; 

I=8.46825e-12; 

E=71e9; 

n=8; 

dampf =0.0036; 
modes ) 

f=zeros ( 2*n, 1 ) ; 
f(2)=l; 

kh=0. 182926; 

T=0 .01; 


%  linear  mass  density  (kg/m) 

%  lumped  mass  (kg) 

%  lumped  mass  rotary  moment  of  inertia  (kgm  2) 

%  elemental  length  (m) 

%  1st  element  length  (m) 

%  8th  element  length  (m)  ^  .  a.  •  /i,  -on 

%  beam  cross  sectional  moment  of  inetria  (Kgm  ) 

%  Youngs  modulus  (aluminum) 

%  #  of  elements  o 

%  Estimated  damping  ratio  for  beam  (Exp  from  1st  2 

%  force  vector  (actuator) 

%  H  infinity  gain  (due  to  beam  properties) 

%  Sampling  period  for  discretized  system 


%  Elemental  mass  and  stiffness 

mll=m*h*[156  22*h;22*h  4*h''2]/420; 
m22=m*h*[156  -22*h;-22*h  4*h“2]/420; 
ml2=m*h*[54  -13*h;13*h  -3*h“2]/420; 
m21=ml2 ' ; 

mlle8=m*h8*[156  22*h8;22*h8  4*h8'2]/420; 
m22e8=m*h8*[156  -22*h8 ; -22*h8  4*h8  2]/420; 
ml2e8=m*h8* [54  -13*h8;13*h8  -3*h8  2]/420; 
m21e8=ml2e8  * ; 

%  Global  mass  matrix 

gm(l:2,l:2)=m22+mll; 
gm( 1 : 2 , 3 : 4 ) =ml2 ; 
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for  i=l:n-2 

gin(2*i+l :  2*i+2, 2*i-l :  2*i)=m21; 
gin(2*i+l:2*i+2,2*i+l:2*i+2)=m22+mll; 
gin(2*i+l:  2*i+2, 2*i+3  :  2*i+4)=inl2; 

end 

gin(2*n-l:2*n,2*n-3:2*n-2)=m21; 
gin(2*n-l :  2*n,  2*n-l :  2*n)=in22; 


gin(13:14, 13:14)  =m22+inlle8 ; 
gin(  15: 16, 13: 14)=m21e8; 
gIn(13:14,15:16)=ml2e8; 
gIn(15:16,15:16)  =m22e8 ; 

%  Add  lumped  mass  (possible  9  nodes) 

MI=[ml  0;  0  II]; 

node=8 ; 

p=2* (node-1 ) ; 

gm(p-l:p,p-l:p)=gin(p-l:p,p-l:p)+MI; 

%  Elemental  stiffness  matrices 

kll=E*I*[12  6*h;6*h  A*h'2]/h^3; 
k22=E*I*[12  -6*h;-6*h  4*h“2]/h~3; 
kl2=E*I*[-12  6*h;-6*h  2*h''2]/h''3; 
k21=kl2'  ; 

klle8=E*I*  [12  6*h8;6*h8  4*h8''2]/h8'‘3; 
k22e8=E*I* [12  -6*h8;-6*h8  4*h8 ^ 2 ] /h8 “ 3 ; 
kl2e8=E*I* [-12  6*h8;-6*h8  2*h8“2] /h8''3; 
k21e8=kl2e8' ; 

%  global  stiffness  matrix 

gk( 1:2, l:2)=k22+kll; 
gk( l:2,3:4)=kl2; 

for  i=l:n-2 

gk(2*i+l:2*i+2,2*i-l:2*i)=k21; 

gk(2*i+l:2*i+2,2*i+l:2*i+2)=k22+kll; 

gk(2*i+l:2*i+2,2*i+3:2*i+4)=kl2; 

end 

gk(2*n-l:2*n,2*n-3:2*n-2)=k21; 
gk(2*n-l: 2*n, 2*n-l: 2*n)=k22; 


gk( 13: 14, 13:14)=k22+klle8; 
gk( 15: 16, 13: 14)=k21e8; 
gk(13:14,15:16)=kl2e8; 
gk(15:16,15:16)=k22e8; 
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%  Find  mode  shapes  and  natural  frequencies 


[V,D]=eig(gk,gm) ; 

[ lambda, k]=sort(diag(D) ) ; 
V=V(:,k); 


factor=diag(V' *gm*V) ; 
phi=V*inv(sqrt(diag(factor) ) ) ; 
omega=diag( sqrt (phi ' *gk*phi ) ) ; 
f req=omega/2/pi ; 
modf=phi ' *f ; 


%  Calculate  the  open  loop  transfer  function  of  the  beam  at 
%  each  frequency  for  the  first  md  modes  (truncate). 

md=input (' Number  of  modes  (1<#<8):  '); 

s2=diag(ones (md) ) ; 
sl=2*dampf *omega( 1 :md) ; 
sO=omega ( 1 : md ) . *omega ( 1 : md ) ; 

num=phi(2, 1 :md) ' . *modf ( 1 :md) ;  %  sensor  measures  theta  #  2 

den=[s2  si  sO]; 

[numg,deng] =parallel (num( 1 ) , den( 1,1:3), num( 2 ) , den( 2 ,1:3)); 
for  i=3:md 

[numg,deng]=parallel(numg,deng,num(i) ,den(i, 1:3) ) ; 
end 

%  Discretize  open  loop  system 
[numgd,dengd]=c2dm(numg,deng,T, ' tustin' ) ; 


%  Control  Transfer  Function 

flagl='y ' ; 
while  flagl=='y' 

disp('(l)  Integral') 

disp('(2)  Derivative  &  Filter') 

disp ( ' ( 3 )  sqrt ( s )  '  ) 

disp('(4)  sqrt(s)  &  Filter:  ') 

case=input (' Enter  Controller  Type:  ' ) ; 

gain=input (' Enter  gain  factor:  '); 

%Estimate  square  root  of  s 
zrs  =  [-.0001  -.01  -1  -100]; 
pis  =  [-.001  -.1  -10  -1000]; 
nums=poly ( zrs ) ; 
dens=poly (pis ) ; 
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[nuinsd,densd]=c2dni(nuins,dens,T,  'matched'  )  ; 

%  Noise  filter 

wc=input( 'Enter  Filter  cutoff  frequency  (rad/sec):  ); 

dmpf=. 707; 

numf=wc"2; 

denf=[l  2*dmpf*wc  wc''2]; 

[numfd, denfd] =c2dm( numf , denf / T, ' tustin ' ) ; 

if  case  ==  1 

%  Integral  control 
numhd=gain* [ T  0 ] ; 
denhd=[l  -1]; 
elseif  case  ==  2 
%  Derivative 

numhd=conv ( numf d , [ 1  - 1 ] ) ; 
denhd=conv( denfd, [T  0]); 
elseif  case  ==  3 
numhd=gain*numsd ; 
denhd=densd; 

elseif  case  >=  4  _  .  j  jv 

[ numhd , denhd ] =ser ies ( numf d , denfd , gain*numsd , densd ) ; 

end 

%  Closed  loop  system 

[ numc Id , dencld ] =f eedback ( numgd , dengd , numhd , denhd , -l ) ; 


%  Frequency  Domain  Evalutaion. 

%  Plot  Bode  for  open  and  closed  loop  systems. 

w=logspace ( 0 , loglO ( 600 ) , 200 ) ; 

[ mg , pg ] =dbode ( numgd , dengd ,T,w); 

[ me , pc ] =dbode ( numc Id , dencld , T , w ) ; 

figure(l) 

semilogx(w, 20*logl0 (mg ) , w, 20*logl0 (me ) ) 

title (' Discrete  Bode  Plot:  Open  and  Closed  Loop  (Ts-O.Ol  sec)  ) 
xlabel ( ' omega  ( rad/sec ) ' ) 
ylabel( ' |G(s) | ' ) 

figure ( 2 ) 

semilogx ( w , pg , w , pc ) 

title (' Discrete  Bode  Plot:  Open  and  Closed  Loop  (Ts=0.01  sec)') 
xlabel ( 'omega  (rad/sec) ' ) 
ylabel( 'Phase' ) 


%  Time  Domain  Evaluation 
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%  Plot  output  of  open  and  closed  loop  system. 


flag='y ' ; 
while  flag=='y' 

dist=input( 'Enter  disturbance  frequency  (Hz):  '); 

tmax=20/dist; 

nt=0:T:tmax; 

u=0.05*sin(dist.*6.28*nt) ; 

[yg,xg]=dlsim(numgd,dengd,u) ; 

[ycl , xcl ] =dlsim(numcld, dencld, u) ; 

figure(3) 

plot(nt,yg,nt,ycl) 

grid 

title ( 'Response  of  System  to  Sinusoidal  Input  (Ts=0.01  sec)') 
xlabel('time  (seconds)'); 
ylabel (' theta  (rad)'); 

flag=input ( 'Try  another  disturbance  frequency? (y/n)  ','s'); 

end 

flagl=input( 'Try  another  control  law  /  filter? (y/n)  ','s'); 

end 
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APPENDIX  B 


Finite  Element  Model 


Natural  Frequencies 


l.Oe+003  ♦ 

0.00101324976094 

0.01775301423174 

0.04638342775662 

0.07923661265561 

0.11935988578432 

0.22849778748830 

0.37813429140721 

0.56949274733975 

0.61770877977266 

0.79596117358377 

1.16976439767046 

1.54850209877084 

2.04634262642825 

2.68166333705337 

3.45392614581641 

4.22726231062265 


» 


1  0e+003  * 

Columns  I  through  7 


.006 


(pej)d|Suv 


-.006 


s»|OA 


H-infinity  control  using  200  Hz  sampling  rate 


-.001 


(pej)9|6uv 


siiOA 
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